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Abstract 

Issues relating to extensive air showers observation by a space-borne 
fluorescence detector and the effects of clouds on the observations are 
investigated using Monte Carlo simulation. The simulations assume the 
presence of clouds with varying altitudes and optical depths. Simulated 
events are reconstructed assuming a cloud-free atmosphere. While it is 
anticipated that auxiliary instruments, such as LIDAR (Light Detection 
And Ranging), will be employed to measure the atmospheric conditions 
during actual observation, it is still possible that these instruments may 
fail to recognize the presence of a cloud in a particular shower observation. 
The purpose of this study is to investigate the effects on the reconstructed 
shower parameters in such cases. Reconstruction results are shown for 
both monocular and stereo detectors and for the two limiting cases of 
optically thin, and optically thick clouds. 



1 Introduction 



Space-borne cosmic rays detectors for energies E > 10 20 eV have been proposed 
[1] and are now under study [2]. Such a detector will comprise one or two 
satellites orbiting the Earth at an altitude of ~ 400 km to ~ 1000 km and will 
have a wide field of view (FOV), on the order of 60°. The footprint on the 
Earth's surface of the FOV has dimensions on the same order of magnitude as 
the orbit height. Studies of the global distribution of clouds and their frequency 
of occurrence, e.g. [3], suggest that the target volume will at any point in time 
contain some clouds. 

The amount of clouds (fractional cover) , the distribution of clouds in terms 
of cloud type, altitude, and optical depth will undoubtedly affect the detector's 
trigger aperture. In addition, cloud presence could result in a reduction of the 
reconstructible aperture, as contaminated events are excluded from the analysis. 
Finally cloud presence could compromise the accuracy of the energy estimate 
for an observed event, since this estimate depends in part on a knowledge of the 
atmospheric conditions at the time and location of the shower development and 
along the path the light from the shower travels to the detector. 

The effect of cloud presence on the detector aperture is beyond the scope of 
this paper. In this study we limit our attention to the question of how cloud 
presence may affect the reconstructed shower geometry and energy. In the con- 
text of a Monte Carlo study, this question can be addressed by applying the 
event analysis assuming no cloud presence, and then determine (a) whether 
the reconstruction procedure can identify the presence of otherwise unreported 
clouds, and thereby rejecting the event in question, and (b) for those events 
where clouds eluded all detection attempts, how the reconstructed shower pa- 
rameters were altered. 

With respect to the detector itself, there are two possible modes of operation: 
monocular and stereo, the latter employing two sites (satellites) separated by 
some distance and which view the same region of the sky. The Fly's Eye experi- 
ment has demonstrated the superiority of the stereo technique on the ground [4] . 
For space-borne detectors, it has been suggested [5] that monocular observation 
can perform as well as stereo if use is made of the information provided by the 
reflection of the Cerenkov beam associated with the shower off the surface of 
the Earth, in order to reconstruct the shower geometry. In this study we also 
investigate possible errors introduced in cases where the reflection occurs off the 
top of a cloud instead of the surface. 

The answer to (a) above will depend on whether or not cloud presence will 
manifest itself through a significant alteration in the expected detector response 
to the shower signal. As an example, the reflectivity of an optically opaque 
cloud is several times larger than that of the surface of the ocean, ~ 80%- 
90% vs. ~ 10%-20%, therefore a test may be be developed and applied to an 
individual shower observation looking at the signal strength of the last few pixels 
to determine whether the reflection of the beam has occurred off the top of a 
cloud. The development of such a test is not trivial. It must be applicable to a 
wide range of shower energies and geometries as well as accommodate different 
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atmospheric conditions and cloud optical properties. Also, the formulation of 
such a test must rely on a detailed description of the event data recorded by 
the detector for each shower observation. As described in section 3.1, we do not 
attempt a detailed simulation of the detector data acquisition system and event 
formation logic. Also, we only treat a few combinations of shower geometries and 
clouds configurations. Therefore, the development of a test for cloud presence 
based on the event data is beyond the scope of this study. 

Clouds come in a wide variety of cloud types, heights, vertical extent, and 
optical depths. There are, in general, also spatially in-homogencous and finite, 
a few kilometers in lateral extent, clouds. This makes a general treatment of 
all possible scenarios difficult. To simplify the discussion we will concentrate on 
two limiting cases. The first case is that of a high altitude, optically thin cloud. 
This case corresponds to cirrus clouds which are pervasive in the atmosphere 
[3]. The second case is that of low altitude, optically thick clouds. These types 
of clouds are easy to detect in general but may be difficult to detect under some 
circumstances, e.g., if the cloud is small in lateral size (on the order of a few 
kilometers.) 

This paper is organized as follows: The next section provides a motivation 
for the different cloud configurations used in the study. Following that is an 
overview of the Monte Carlo simulation of the detector, showers, and the atmo- 
sphere including cloud simulation. Section 4 provides an overview of the shower 
geometry and energy reconstruction procedures. Finally section 5 presents the 
results of the study. 

2 Clouds and EAS 

Clouds are classified as (a) low-level (cloud base height, hb ase < 2 km), (b) mid- 
level (2 < hb ase < 6 km), and (c) high-level ( hbase > 6 km) [6]. In the equatorial 
region, high-level clouds (Cirrus) typically occur at altitudes of 8 — 15 km [6]. 
Most EAS develop in the lower atmosphere at altitudes < 20 km, where most of 
the atmospheric mass is located. Depending on their altitude, clouds are made 
up of predominantly water molecules, (low- level), mixed water molecules and 
ice crystals (mid- level), and ice crystals (high-level). For our purposes, water 
and ice crystal clouds have one important difference: the number density of 
scatterers. The concentration of water molecules in a low level cloud is a factor 
of 10 to 100 greater than that of ice crystals in a high altitude cirrus cloud. 
This results in a much smaller scattering length for the low level cloud, and 
the relation between the optical depth of the cloud and its physical thickness 
becomes qualitatively different. 

From detailed Monte Carlo simulations, most EAS generated by protons or 
nuclei in the energy range 10 19 < E < 10 21 eV reach maximum development 
at atmospheric slant-depths, x, between 700-1000 gm/cm 2 (depending on the 
energy, primary type, and the hadronic model used in the simulations) [7]. 
Beyond the shower maximum depth, x max , the number of electrons in the shower 
falls rapidly. 
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Table 1: The height (km) above the surface, of a point along the shower track 
at a given slant depth (gm/cm 2 ) along the shower, for showers with different 
zenith angles, 0. Atmospheric density profile according to the U.S. Standard 
Atmosphere, 1976. 
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To quantify how clouds might affect space-borne observations of extensive 
air showers, we need to relate the atmospheric slant depth along the shower 
track to altitude above the Earth's surface (sea level). This relation depends on 
the shower zenith angle and is presented in figure 1 and table 2. As can be seen 
from the table, different cloud altitude and shower zenith angle combinations 
can result in the shower front reaching the cloud top at different stages of the 
shower development. Three broad cases can be identified: (1) clouds above 
shower development, (2) clouds in the region of shower development, and (3) 
clouds below shower development. 

In addition to the location of the cloud in relation to the shower-detector 
geometry, the optical depth, r, and physical thickness, Az, of the cloud also 
play a role in determining what the detector sees. Hence they also need to 
be considered in combination with the height of the cloud top, h top . Finally, 
depending on the values of r and Az, the effects of multiple light scattering may 
or may not be negligible, <10%. 

As mentioned in the introduction, we will only address the two limiting cases 
for the optical depth. The reason for this is that the consideration of these cases 
is sufficient for the study of the relevant problems of: light transmission through 
high altitude cirrus clouds (optically thin) , and the Cerenkov beam reflection off 
optically thick clouds and its effect on the monocular geometry reconstruction. 
In general, single scattering calculations are sufficient for a medium r < 0.1. 
First order corrections may be required for 0.1 < r < 0.3 [8, 9]. From our own 
studies, we saw that for r < 0.5, second order corrections account for less than 
~ 10% of the total signal transmitted through a cloud. Therefore, to avoid 
having to calculate light multiple scattering beyond first order corrections we 
will restrict our definition of optically-thin clouds to mean clouds with optical 
depth t < 0.5. 

For the case of an opaque cloud, a light beam impinging on the top of the 
cloud will be reflected as a result of a large number of multiple scatterings inside 
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Figure 1: Slant depth (g/cm 2 ) along a shower track with a zenith angle, 6, vs. 
the height (km) above the surface of a point along the shower track. Data from 
table 1. 

the cloud. The amount, spatial and temporal distributions, and the direction 
of the reflected photons can be calculated using a Monte Carlo procedure (see 
section A. 3). The use of this procedure or an equivalent detailed simulation of 
the cloud reflective properties would be required if one is to attempt to infer 
the cloud presence from the event data. For this study, however, it is sufficient 
to treat reflection off clouds in the same fashion as reflection off the surface, as 
described in section 3.4. 

3 Simulation 

Both the simulation and reconstruction programs used in this study arc based 
on those developed for the High Resolution Fly's Eye (HiRes) experiment [10]. 
The original programs are described in detail in [11]. While the underlying 
algorithms arc similar to those used by HiRes, the actual code was converted 
from Fortran and C to the C++ language and extensive use was made of the 
ROOT data analysis framework [12]. 

Naturally, the HiRes detector simulation was replaced by a description of 
the OWL detector. Otherwise, a large portion of the simulation code, e.g. 
the atmosphere, is detector independent, and was retained. Cloud simulation 
including the effects of multiple light scattering was added, and minor modifi- 
cations were made throughout the code to account for the differences between 
the two detectors. 

In the following subsections a description of the detector simulation is given, 
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followed by a synopsis of the atmospheric modeling. Wc then present an overview 
of the shower simulation, and the section ends with a discussion of the Cerenkov 
spot. 

3.1 The Detector 

A description of the OWL baseline instrument is given in [13]. Earlier versions 
of proposed designs were presented in a workshop [2] . Our detector simulation 
is based on these earlier designs. The conclusions drawn from these studies 
will not be substantially affected by the design evolution. A description of the 
simulated detector follows: 

The detector consists of two cameras, each mounted on a satellite. The 
satellites have an orbital height of 800 km and are separated by a distance of 
500 km. Each camera comprises a concave, spherical mirror and a focal plane 
detector. The mirror has an effective light collecting area of 4.9 m 2 and a field 
of view (FOV) corresponding to a cone of half-angle of 30°. The mirrors axes 
are tilted slightly from the nadir in order for the two cameras to view a common 
area on the surface. 

For the reason that the optical design of the detector had not been completed 
at the time this study was begun, and also for the sake of simplicity, we opted 
to use a scheme in which all photon ray-tracing calculations are done in angular 
space. In this case we ignore the details of the detector optics and simply treat 
the mirror as a "light collector" with a circular aperture. The focal plane pixels 
are arranged on a rectangular grid with each pixel having a fixed angular size. 
The pixel angular size is selected to meet the detector design requirement of 
resolving a distance of 1 km on the surface. So, for an orbit height h, the pixel 
angular size, <5, is equal to 1/h radian, with h measured in km's. For a 800 km 
orbit this translates to 5 w 0.7°. We assume full coverage of the focal surface, 
i.e., we ignore the physical gaps and dead areas between the PMT's. 

During an event simulation we calculate and record the arrival time of each 
photon reaching the detector from the shower. The arrival direction of the 
photon determines the pixel in which the photon is registered. Data for a 
triggered pixel comprises the pixel pointing direction, the integrated pixel signal 
and the mean arrival time of the photons recorded by the pixel. The integrated 
signal is simply the number of photo-clcctrons (pe) recorded by the pixel in a 
time window of 12 /is. The particular choice for the width of the time window 
allows enough time for a shower with a zenith angle of 30° or greater to cross 
the field of view of the pixel. 

Wc employ a simple detector trigger scheme, which requires at least six pixels 
to fire from the light of the shower in order to form an event. A pixel trigger 
occurs if the pixel records three or more pe in a one fj,s interval. The test for an 
individual pixel trigger is performed as follows: 

1 . the arrival times for each of the pe recorded by the pixel are sorted in time 
to find the arrival times of the first and last recorded pe A time window 
is formed around these times, and a 2 fis interval is added to each end. 
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2. sky noise pe are added to the expanded time window assuming a uniform 
background of 200 photons/m 2 /sr/ns. 



3. shower generated and noise pe are now stored in a histogram with a 1 /is 
bin width. The histogram's bins are scanned and if any bin has three or 
more entries then the pixel trigger flag is set. 

4. finally, if the total width of the time window exceeds 12 /is, then a sliding 
window of that width is used to scan the histogram to find the set of 
contiguous bins with the largest sum. 

Finally, we incorporate elements from the HiRcs detector to cover some of the 
gaps in the simulation. In particular we assume that the detector will use a UV 
filter similar to that of HiRes, which passes light in the 300-400 nm range. A 
parameterization of the wavelength dependence of the PMT quantum efficiency 
is also borrowed from the HiRcs simulation, so is a constant mirror reflectivity 
of 80%. These detector components are described in detail in [10]. 

3.2 The Atmosphere 

There are four elements or components to the simulation of UV light transmis- 
sion through the atmosphere. These include Rayleigh scattering by molecules of 
air, scattering by surface aerosols, absorption by ozone molecules, and scattering 
by clouds. The treatment of the first three is based on the HiRes simulations; 
we adopted the same models without modifications. Although these models are 
more appropriate for the Utah desert observation conditions than for observa- 
tion over the ocean, the differences should have little effect on the results of 
this study. This is because the shower development and light propagation to 
the satellites occur almost entirely above the surface aerosol layer, which is the 
one factor most likely to be significantly different between the desert and the 
ocean. Before we turn to a discussion of the cloud simulation we present a brief 
overview of these models. 

Light scattering is characterized by the scattering cross section, (5, and the 
phase function, P = P(cos8 s ), where 9 S is the scattering angle. The cross 
section for molecular, or Rayleigh, scattering is given by: 



where /3r is measured in units of m _1 , p(h) is the air density (g/cm 2 ) at altitude 
h (m) above sea level, and xr = 2970 g/cm 2 is the mean free path at wavelength 
A = 400 nm. The air density and temperature profiles as function of altitude 
are given by the U.S. Standard Atmosphere, 1976 [14]. The phase function for 
Rayleigh scattering is given by: 



fa = 100( /9 (/ l )/a; i? )(400/A) 4 



(1) 




(2) 
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Aerosols scattering is calculated according to the following formula: 




(3) 



Here L a is the scattering length at the surface, and we have: 



/?„ (A) = Pa /L a {\) 
The aerosols reduced density, p a , is given by: 



(4) 




1 h < h 

e -{h-h m )/H a h > h 



(5) 



where h m is the height of the mixing layer, and H a is the scale height above 
the mixing layer [15]. The scattering length at wavelength A = 334 nm is 
a free parameter of the model. The wavelength dependence of the scattering 
process is accounted for by a parameterization, shown in figure 2. The aerosols 
scattering phase function, also shown in figure 2, is based on the Longtin desert 
aerosols model [16]. For this study, the model parameters are set to: h m = 0, 
H a = 1.2 km, and L Q (334) = 23.0 km. 

Ozone absorption is characterized by a wavelength-dependent absorption 
coefficient, 0*03, and an altitude-dependent concentration, po3, both shown in 
figure 2. Model parameters were extracted from [17]. The extinction length, in 
meters, due to ozone absorption can be written as: 



with the constant factor accounting for unit conversion. 

The cloud model and the simulation of light propagation in a cloud is de- 
scribed in detail in appendix A. Here we present a brief description of the 
model. 

Simulated clouds have a uniform density which steps to zero at the cloud 
boundaries. The parameters used to describe a cloud are the cloud base height, 
htase, cloud top height, h top , and the optical depth r. The scattering length, {3, 
inside the cloud is related to the optical depth by the relation: 



The wavelength dependence of the scattering process, which is mild in the near 
UV [18], is ignored. Finally, scattering within optically thin clouds calculatioins 
use a phase function appropriate for ice clouds at A = 0.5 /im [19]. 

3.3 Shower Simulation 

The primary particle energy is selected at the start of the simulation. The shower 
track geometry is generated randomly in order to obtain uniform and isotropic 
showers distributions in the atmospheric volume viewed by the detector. 
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(9.87 x 10~ 7 ) x a 03 (A) x P03 (h) 



(6) 



L03 



T = P{htop - hbase) 



(7) 
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Figure 2: From the top, the wavelength dependence of the aerosols extinc- 
tion length, the aerosols phase function, the wavelength dependence of ozone 
absorption, and the ozone concentration as a function of altitude. 
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The generated shower profiles follow a Gaisser-Hillas function [20]. The 
number of electrons, N ei in the shower is given as a function of depth, x, by: 



N e (x) = N max X ~ X _° e (x mM -x)/A (8) 

where the Xq parameter is chosen from an exponential (< xq >=35 g/cm 2 for 
proton primaries), x max is chosen from a Gaussian distribution, A is fixed at 70 
gm/cm 2 , and N max is selected so that the integral of the profile, corrected for 
lost energy [21] [22], gives a total shower energy equal to that of the primary 
particle: 

Etot = E corr + 2.18 x J N e dx (9) 

The mean and variance of the x max Gaussian distribution depend on the 
energy and mass number of the primary cosmic ray particle. For protons we 
assume < x max >= 725 g/cm 2 at E = 10 18 cV, increasing by 55 g/cm 2 per 
decade in energy. The standard deviation is set to 50 g/cm 2 for all energies. 
These values are based on shower simulations quoted by the Fly's Eye group in 
their analysis [23] , and are consistent with simulations results from the CORSIKA 
program with the QGS Jet model [7] . 

The Nishimura-Kamata-Greisen (NKG) function [24, 25], is used to describe 
the lateral distribution of shower electrons: 

(\ N f( r 



where 



s-2 / \ s-4.5 

r 



f _ r(4.5 - s) f 

27rT(s)T(A.5-2s)\r M J 

and the shower age parameter, s, is given by: 



(10) 



s = 3x/(x + 2x max ) (11) 

Here, tm is the Moliere radius. The value of tm depends on the air density and 
is evaluated at each point along the track. 

Fluorescence light is generated according to the formulas given in [15], but 
more recent measurements of the air fluorescence yield are used [26] . The calcu- 
lation of the Cerenkov light production also follows that of the Fly's Eye paper 
[15]. 

The above procedure is modified if optically thin clouds are present. First 
wc identify all track segments which lie inside the cloud. For these segments 
the simulation proceeds as described above but in addition, the number of fluo- 
rescence photons scattering once in the cloud, and Gerenkov photons scattering 
once or twice in the cloud, before reaching the detector, are calculated. These 
photons are included in the detector response as additional signal. 
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3.4 Cerenkov Spot 

The Cerenkov spot refers to the area on the surface around the shower core 
where the Cerenkov beam is reflected, off the surface and into the detector. 
The simulation of the signal recorded by the detector, and generated by the 
reflected beam, requires the consideration of three separate issues: (1) The 
lateral distribution of the Cerenkov photons at a given point along the shower 
development, (2) the lateral distribution of the Cerenkov photons at the surface, 
and (3) the reflection properties of the surface. 

The lateral distribution of photons in the Cerenkov beam is assumed to fol- 
low that of the shower electrons, i.e. it's given by the NKG function. This is 
a simplification and in general results in a greater concentration of Cerenkov 
photons near the shower axis. A more accurate description of the lateral spread 
of the Cerenkov beam would be required ( along with a detailed detector sim- 
ulation ) to address the problem of identifying cloud presence from event data. 
For this study, however, the use of the NKG function is sufficient. 

The Cerenkov front has a circular shape centered around, and perpendicular 
to the shower axis. The spot formed on the surface by the beam is in general 
elliptical with the elongation of the spot depending on the zenith angle of the 
shower. In the simulation, the transformation of the photon position from a 
point on a circular disk about the shower axis to a point on the reflecting surface 
is performed and the time offset of the photons relative to the shower core is 
calculated before the photon is propagated to the detector. 

Finally, the surface reflection albedo is assumed to be constant at 20% for 
reflection off water, and the reflection is assumed to be isotropic (into 2ir). The 
same is assumed for cloud reflection. 

4 Event Reconstruction 

Events are reconstructed from the raw event data to obtain the shower energy, 
shower x max , and the arrival direction of the primary cosmic ray particle. The 
raw data consists of a set of triggered pixels with known pointing directions, each 
with a measured mean arrival time and a time- integrated total pe count. With 
two instruments observing the shower, the combined data from each instrument 
comprises a stereo event. Data from each instrument can be analyzed separately 
as a monocular event. 

The reconstructed shower is described by a set of parameters which specify 
the shower geometry and shower profile. The shower energy is obtained from 
the shower profile using cq. 9. The shower track geometry can be described by 
a pair of orthogonal vectors R p and Ut, the latter being the shower direction 
unit vector. Alternatively the geometry can be specified in terms of the Shower- 
Detector (SD) plane normal, h, and a pair of scalars R p and ip which determine 
a line in that plane. The shower profile is given by the Gaisser-Hillas function, 
eq. 8. 

Event reconstruction is divided into three consecutive steps: 
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1. SD plane reconstruction for each eye. 

2. Shower track geometry reconstruction. 

3. Shower profile and energy reconstruction. 

Step 2. above is implemented differently for monocular and stereo events. For 
stereo, the intersection of the two SD planes from each eye determines the shower 
track. Monocular reconstruction requires the use of pixel trigger timing and an 
additional constraint provided by the observation of the Cerenkov spot. All 
other steps are similar for both monocular and stereo events. 

In general, a file containing a set of Monte Carlo generated events contains 
a reference to the set of atmospheric parameters used in the simulation. This 
enables the reconstruction programs to use the same atmosphere used in the 
simulation. However, since the purpose of this study is to investigate the effect 
of clouds which go undetected, on the event reconstruction, clouds are removed 
from the atmosphere during reconstruction. 

4.1 Shower Geometry 
4.1.1 Shower-Detector Plane 

The shower-detector plane is that plane which contains the detector, a point, 
and the shower track, a line. For a stereo detector a SD plane is calculated 
separately for each eye. The SD plane is calculated by minimizing a \ 2 function 
given by: 



where the sum is over triggered pixels, n is the plane normal, hi are the pixel 
viewing direction vectors and Wi are weights equal to the total number of pc 
seen by pixel i. An angular pointing error of a ~ 0.07° (equal to the pixel 
angular size) is assumed for all pixels. 

4.1.2 Shower Track in the SD plane 

In the case of stereo observation, the intersection of the SD plane normals from 
each eye describes a line in space, namely the shower track. This method despite 
its simplicity works very well in general [27]. Only events for which the opening 
angle between the two planes is small and the plane determination was not good, 
e.g. due to short track-length, does the method fail to produce accurate results. 

In the case of monocular observation, track reconstruction uses the pixel 
timing information. The timing fit method is based on the relation between the 
crossing time of the shower front in a pixel's field of view (mean photon arrival 




(12) 
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time at the detector) and the pixel's viewing angle. The pixel crossing time U 
as a function of the pixel viewing angle Xi m the SD plane, is given by: 

R 1 

U = t H ?-tan-(ir - - Xi) (13) 

The reader can refer to [15] for a derivation and definition of the parameters. 
In all, there are three unknown fit parameters, namely t Q , R p , and ip. 

A fit based on eq. 13 produces accurate results only when the range of an- 
gles covered by the triggered pixels, i.e. the angular track-length of the event, is 
"large enough". A good discussion of the timing relation and the requirements 
for accurate reconstruction appear in [28] in relation to the Auger detector. For 
showers observed from space, the angular track-length is too small for the fit to 
result in satisfactory results, and an additional constraint on the shower geom- 
etry is required. The observation of the Cerenkov spot provides this constraint 
in the form of a known shower impact point on the surface, also referred to 
as the shower core position. The use of this constraint results in a significant 
improvement in the accuracy of the fit. The shower core vector is obtained from 
the event data as follows: 

If the reflected Cerenkov light is observed by one pixel with a pointing direc- 
tion vector v c then, the shower core vector can be calculated using the relation 
r c = r m + sv c where r c is the shower core position, f m is the detector position 
vector, and s is a scalar which can be solved for by making the requirement 
that this line intersect the surface of the Earth (a sphere with a known radius 
corresponding to an altitude of h = m above sea level). 

In general, the reflected light is observed by one or more pixels, to identify 
which we examine the set of triggered pixels for the one that triggered last in 
time and the one triggered last in angular distance from the start of the track. A 
group of one or more pixels is first identified as being triggered by the reflected 
beam by examining those pixels adjacent to the last triggered pixel for their 
trigger times. A sum of the pointing directions of these pixels (weighted by 
the total signal in each pixel) is performed to get an average direction. This 
direction is then projected in the SD plane to get a final estimate of the core 
direction. 

4.2 Shower Profile 

With the shower track geometry in hand we proceed to reconstruct the shower 
profile. The shower profile is assumed to follow a Gaisser-Hillas function with 
three free parameters: a;o, x maxi and N max . For each trial profile a shower 
is generated with the reconstructed geometry and the detector response to the 
shower is calculated. The calculation proceeds along the same lines as the Monte 
Carlo (same light production and propagation models) with the exception that 
all random fluctuations are suppressed. Where in the MC the number and start- 
ing position of ray-traced photons is chosen randomly; During reconstruction, 
the mean number of photons from each track segment is distributed on a two 
dimensional grid representing an NKG lateral distribution. The effect of the 
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finite mirror spot size is accounted for by distributing the flux received by the 
detector among the pixels according to the distribution of the spot. 

The best fit (reconstructed shower profile parameters) is chosen to minimize 
a x 2 function calculated from the observed (MC output) and fit values for the 
pe counts from each pixel in the event. The following function is used: 

i * 

where the sum is over triggered pixels, is the measured pixel signal in pe, 
is the predicted pixel signal, and of = + £?,. The of terms are obtained 
by adding in quadrature the Poisson fluctuation in the signal, vSH, and the 
estimated sky background fluctuations for that pixel. It should be noted that 
not all triggered pixels are included in the sum. Those pixels near the end of 
the track believed to be triggered by the reflected Cerenkov beam are excluded 
from the sum. 

5 Results and Conclusions 

In this section we summarize the results obtained from the two studies of the 
possible effects of cloud presence on the reconstructed shower energies and x ma x ■ 
We start with optically thin clouds. 

5.1 Optically Thin Clouds 

The effect of cloud presence is examined as follows: A Monte Carlo shower is 
generated in a cloud free atmosphere and the detector response is evaluated. 
Next, a loop over a set of cloud configurations (described below) is made in 
which the selected cloud is included in the atmosphere simulation. The same 
shower from the cloud free simulation is developed through the atmosphere, and 
the detector response is recorded. All events are reconstructed assuming a cloud 
free atmosphere. 

Of the many (infinite) possible cloud configurations we selected the following 
set: Cloud base height is set to 6 km for all clouds in the study. Four different 
cloud top heights are used, these are given in table 2 along with the correspond- 
ing atmospheric slant depths along shower tracks at 60° and 75°. At each cloud 
top height setting, the cloud optical depth is varied between 0.1-0.5 in steps of 
0.1, for a total of 20 cloud settings. Most Proton initiated showers with an en- 
ergy of 10 21 eV are expected to have x max values in the range of 800-1000 g/cm 2 . 
If these showers develop at a zenith angle of 60° then they will reach maximum 
development at an altitude just above the selected cloud base height. In most 
cases then the shower will traverse the cloud while it is still increasing in size. 
At 9 = 75°, A cloud top height of 7.28 km or 9.25 km insures that the cloud lies 
below the shower x max . The other cases cover a larger portion of the shower 
development curve. 
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Table 2: slant depth (gm/cm 2 ) vs. height (km) for showers with zenith angles 
of 60°and 75°. The values shown in the first column are chosen 
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Figure 3: Energy and X ma x shift in the presence of optically thin clouds ( r 
plotted along the x-axis ) for four different cloud top height settings, indicated 
by "ht 200" through "ht 800" . Stereo reconstruction of showers with 9 = 60°. 
The function, exp(— r) — 1, is superimposed on the energy plot. 



Figures 3 through 6 summarize the results for stereo and monocular recon- 
struction of shower energy and x max . In the figures, the cloud optical depth, 
t, is plotted along the x-axis, with r = corresponding to the cloud free at- 
mosphere. Four points for each r > correspond to the different h top values. 
Cloud h top values of 7.28 km and 16.23 km are explicitly indicated on the plots 
by their corresponding slant depth values at 6 = 60°. Each point in the plots 
represents the mean shift and standard deviation for a set of 200 reconstructed 
events. 

In the case of stereo reconstruction, the reconstructed shower geometry is 
unaffected by cloud presence. The cloud affects the amount of light reaching 
the detector from different parts of the shower depending on its position and 
extent. In a couple of cases the effect can be easily understood. For example, 
the case of ht op = 16.23 km (ht 200 in the plot) and shower 9 = 60°, the 
mean reconstructed energy is shifted down by a factor of exp(— r). Another 
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Figure 4: Energy and x max shift in the presence of optically thin clouds ( r 
plotted along the a;-axis ) for four different cloud top height settings, indicated 
by "ht 200" through "ht 800". Stereo reconstruction of showers with 6 = 75°. 
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Figure 5: Energy and x max shift in the presence of optically thin clouds ( r 
plotted along the x-axis ) for four different cloud top height settings, indicated 
by "ht 200" through "ht 800". Monocular reconstruction of showers with 8 = 60°. 
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Figure 6: Energy and x max shift in the presence of optically thin clouds ( r 
plotted along the a;-axis ) for four different cloud top height settings, indicated 
by "ht 200" through "ht 800". Monocular reconstruction of showers with 6 = 75°. 



example is provided by the 9 = 75° showers when the cloud lies well below 
shower maximum (htop= 7.28 km or 9.25 km), in these cases we see very little 
effect on the reconstructed shower parameters. 

The dependence of monocular shower geometry reconstruction on the obser- 
vation and correct identification of the Cerenkov spot makes the interpretation 
of the results more complicated. In general, if the shower geometry is recon- 
structed correctly, then the effect of the cloud on the reconstructed energy and 
Xmax will be similar to the effect on stereo reconstructed events. Otherwise, 
the error will depend mostly on what the reconstructed geometry turns out to 
be. Looking at the 8 = 60° showers, one can see that the energy plot looks 
almost identical to that of the stereo case, and that the mean shift in x max is 
comparable but smaller than that for the stereo case, however, the error bars 
arc slightly larger in the monocular case. For most showers in this group, the 
detector did trigger on the reflected Cerenkov light and the geometry recon- 
struction procedure gave the right results. In a few cases, as the cloud optical 
depth increased, the detector did not trigger on the reflected beam and a wrong 
geometry resulted based on the false identification of the ground spot. 

For showers with 6 = 75°, the true Cerenkov spot was not observed in a large 
number of cases. This can be explained by noting that: (a) the larger inclination 
of the shower means that the Cerenkov beam will go through a larger distance 
through the cloud resulting in more attenuation, and (b) especially for the cases 
in which the cloud lies below the shower maximum development, there are no 
more shower particles to feed the Cerenkov beam. Figure 7 shows an example 
of a reconstructed event in the presence of a cloud with r = 0.4 lying below the 
shower development. 
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Figure 7: Example showing the effect of a cloud on the calculated light flux 
(pe/O.ldeg/m 2 ) seen by the detector. In this example a 10 21 eV proton shower 
with a zenith angle of 75° passes through a cloud with an optical depth of 0.4 
which extends between 6 km (w 1500 g/cm 2 ) and 7.28 km (w 1800 g/cm 2 ) 
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Figure 8: Energy and x max shift for 10 21 eV showers with zenith angle 
as a function of cloud top height. 
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5.2 Optically Thick Clouds 

For this study we place optically thick clouds with h top in the range of 1-5 km, 
in one km steps, and generate showers at a fixed energy of 10 21 eV and with 
fixed zenith angles of 60° and 75°. Sets of 500 events each were generated 
for each combination of shower zenith angle, and cloud height. All the data 
sets were reconstructed and analyzed assuming no clouds. A number of quality 
cuts were applied to the reconstructed sets of showers in order to remove badly 
reconstructed events. A list of the applied cuts follows: 

1. xlfi/ndof < 10 

2. SD planes opening angle > 6°. (Stereo events) 

3. Observed angular track-length > 0.6° (approx. 9 pixels) 

4. Number of good angular bins > 5 
5- Xmax bracket cut: Xfi rs t + 100 < 

Xmax • 

Cuts number 2 and 3 remove events where the geometrical reconstruction results 
were probably not accurate. Cut number 4 is somewhat similar to 3, but is over 
good angular bins. The last cut represents the requirement that the shower x max 
was observed. Only events which passed the cuts are included in the results. 

Figure 8 shows the resulting shift in energy and x m ax, for the case 8 = 60°. 
Figure 9 shows the same for 6 = 75°. The results are also shown in tables 
3, and 4. The #events in these tables indicates the number of successfully 
reconstructed events, out of 500. For the set 6 — 75° and clouds at 4 km the 
reconstruction job terminated before the full set was done, and therefore, the 
smaller number of events. 

The results show that while the performance of stereo is stable for different 
cloud heights, monocular reconstruction suffers badly if clouds at altitudes of 
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Energy shift vs. cloud height xmax shift vs. cloud height 
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Figure 9: Energy and x max shift for showers with 9 = 75°and E = 10 21 eV. The 
x-axis is the cloud top height in km's and the y-axis is (E out — Ei n )/Ei n . 
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Table 3: Energy shift (%) versus cloud height for shower energy of 10 21 eV and 
shower zenith angle of 60°. First column is cloud height in km. the following 
two are for stereo while the last two arc for monocular. 
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Table 4: Energy shift (%) versus cloud height for shower energy of 10 21 eV and 
shower zenith angle of 75°. First column is cloud height in km. the following 
two are for stereo while the last two are for monocular. 
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Figure 10: x 2 / n( V distribution from monocular reconstruction of two sets of 
showers: one set with no clouds, the other with clouds at 4 km. Showers were 
generated with a fixed energy of 10 21 eV and a zenith angle of 60°. 



2 km or higher are present and their presence goes unrecognized. The recon- 
structed shower profiles in the monocular case "look" normal with the exception 
of the abnormal development depth and result in reasonable values for the % 2 
as shown, for an example, in figure 10. 

In case of stereo geometry, the last observed point along the shower track 
(pointing direction) can be converted to a position in space and therefore a 
height above the surface. This point can be interpreted as the surface height 
or cloud top height. Cloud presence can be identified by comparing this height 
with the known surface elevation. Figure 11 shows results from a test study. 
The cloud height is underestimated by approximately 0.5 km but the resolution 
is better than 0.5 km. In the figure, the error bars indicate the spread in the 
calculated heights and not the error on the mean. 



21 



Stereo prediction for cloud height 
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Figure 11: Stereo reconstruction of height of shower core, i.e., cloud height 
in case of reflection off opaque cloud. The cloud height is underestimated by 
approximately 0.5 km but the resolution is better than 0.5 km. 
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A Cloud Simulation 

A.l Introduction 

The input parameters for the clouds model are: the cloud base height, hbase, 
the top height, h top , and the vertical optical depth, r. The clouds density is 
uniform and falls to zero at the boundaries: 

1 h base < h < 

'Hop 

otherwise 

The vertical extent of the cloud, Az = h top — h base , and the optical depth, 
r, determine the scattering coefficient through the relation r = /3Az. With (3 
given in 1 /m. Note that a e = (3 since we assume no absorption. The wavelength 
dependence of the cloud's optical parameters is mild for A < 0.5 /im [18] and is 
ignored. 

Several phase functions are relevant to a discussion of cloud scattering. The 
simplest is the Henyey-Greenstein (HG) phase function [29], given in eq. 15. It 
is often used in radiative transfer calculations as an analytic approximation to 
actual phase functions which may display complicated structures, see [30] and 
references therein. 

JWcoB(fl.)!ff) = — (is) 
[1 + g 2 - 2gcos{9 s )\ ' 

Note that the parameter g in the HG function is equal to the asymmetry pa- 
rameter defined by: 

1 f 1 

g = - P(cos(9 s )) cos(9 s )dcos(9 s ) 



where 9 S is the scattering angle and P(cos(0 s )) is the phase function. 

One feature of realistic clouds phase functions is a backward scattering peak. 
This feature is not reproduced by the HG function, however, a double-Henyey- 
Greenstein (DHG) function can provide a better fit, see fig. 12. The DHG 
function is defined by [31]: 

Pdhg(9) = f x P HG (gi) + (1 - /) x Phg(92) 
where / ~ 1 gives the forward scattering strength, and gi is negative. 
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Figure 12: Ice clouds phase function (data from Liou at 0.5 .) A double 
Hcnyey-Greenstein approximation also shown. 

Liou [19] gives in tabular form the phase function for a cirrostratus cloud 
model at a wavelength A = 500 nm. It is shown in figure 12 along with a HG 
and a DHG functions superimposed. A calculation of the asymmetry parameter 
for the realistic phase function gives g = 0.753. The same value is used in 
the superimposed HG function. The DHG function parameters were set to: 
,9i = 0.82, .92 = -0.82 2 , and / = 0.96. 

A. 2 Multiple Scattering in Clouds 

A proper treatment of multiple scattering in clouds has to take into account the 
scattering and absorption by other atmospheric constituents. However, by con- 
sidering the relative strength of the relevant processes, we can show that under 
certain conditions it is safe to ignore some of them. Thin cirrus clouds occur in 
the atmosphere at altitudes greater than 6 km and could be as high as 15 km. 
The extinction length varies with the clouds ice water content (determined in 
part by the altitude) and takes values on the order of a few kilometers. For 
comparison, extinction lengths due to molecular scattering, ozone absorption, 
and aerosols scattering are shown in figure 13. 
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Ozone absorption is negligible for wavelengths greater than A = 320 nm, 
and as can be seen from the figure, it can be safely dropped from the multiple 
scattering calculation for A > 310 nm. Aerosols density and extinction length 
are variable and whether or not they can be neglected depends on the local 
conditions. Figure 13 shows two examples, L a = 10 km, which represents a 
hazy atmosphere, and L a = 23 km, corresponding to an average atmosphere. 
Even in hazy conditions, the attenuation length due to aerosols is large for 
altitudes greater than 6 km because of the small scale height of the aerosols 
density distribution. For an average atmosphere, the aerosols extinction length 
is almost 20 times as large as the typical cloud extinction length, and can be 
safely neglected. Rayleigh scattering has a very strong wavelength dependence 
and can not be ignored at wavelengths close to A = 300 nm. For larger A, and 
at altitudes greater than 6 km, the Rayleigh scattering length is greater than 
20 km. 



A. 2.1 An Isotropic Source 

In this section we develop expressions for the direct transmission and transmis- 
sion due to first, and second order scattering of light from a point source inside 
a cloud to a detector. First we define some notation. The point source is located 
at position fo, the detector (mirror) at f m . The points f\, fi are the locations 
inside the cloud of first and second photon scattering respectively. The distance 
between any of the above points is denoted by l a b where ab take the values of 
the subscripts of the respective points. In the case of the detector, we denote 
by lam the distance between f a and r m which lies inside the cloud, i.e. the 
path-length inside the cloud, the actual distance, \r m — f a \, is denoted by R am - 
Direction and solid angle are denoted by f2. For example, fioi is the direction 
defined by the unit vector (f\ — fo)/ \ f\ — fo\. The detector effective aperture is 
denoted by A, and the projected area of the detector with respect to f by A±\ 
The optical path length due to scattering by non-cloud particles (air, aerosols, 
or ozone absorption) is denoted by r(™ c ) with subscripts to identify the path. 
Non-clouds scattering and absorption coefficients are height dependent and so 
an integration over the path is required. 

Let the point source emit isotropically N s photons, then the number of 
directly transmitted photons reaching the detector is given by: 



iVsA -L „-/3 c i 0m „-T<r> 



First order scattering involves direct transmission to a point r\ then scatter- 
ing in a volume clement dV\ into a solid angle dfl at tt: 



where P(Q, fioi) is the phase function for scattering by the cloud. 
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Figure 13: extinction length due to Rayleigh scattering, ozone absorption, and 
aerosols scattering as a function of altitude above sea level. U.S. 1976 standard 
atmosphere used for density profile for Rayleigh scattering. 
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The number of photons received by the detector due to first order scattering 
taking attenuation along the path from r± to the detector into account, is given 
by the integral over the cloud of: 

The expression for second order scattering is similar to first order, however, 
instead of a point source at r*o we now have an integral over r\ : 



-(r 2 ,fi) = (3 c dV 2 / -^-{ri,n 12 ) 



dn v " ' \] dn 

1 

x p 

'12 



^p{n,n 12 )e~^ h 



The number of photons received by the detector is given by the integral of: 

,(2) 



dN 2 = ^-(f 2 ,n 2m )x ^e-W-e" 
all \ its , 



The evaluation of the above expressions and their integrals is done numeri- 
cally, once a point source, a detector, and the cloud/atmosphere are specified. 

A. 2. 2 Scattering out of a beam 

Instead of an isotropic point source we now consider a light beam propagating in 
the direction, Oo. Here we assume that fio does not point toward the detector, 
i.e. the detector receives no direct light from the beam. Inside the cloud, light 
scattered out of the beam at a point f*o along the beam propagation path and 
into the detector is given by: 

Ni = N s P(n 0m ,n ) x J^e-^e-^ 

\ K 0m J 

where N„ is the number of photons scattered out of the beam at the point ro- 

Second order scattering inside the cloud is treated as follows: First consider 
light scattered out of the beam at ro in some direction f2oi- At a point fi along 
this direction, the irradiance is given by: 

^(fi) = N s p(n 01 ,n ) e -^e- T °" c} 

diiQi 

Next, the scattering in a volume element dV\ at r\ into an arbitrary direction 
n is given by: 



^-(r u n) = p c dv, 
dn 



dN, 

la~ ( ri ) 
diloi 



p(n,n 



01 
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Finally the contribution to the detector signal from the point r\ is: 

The integral over the cloud volume of (IN2 gives the total contribution due to 
second order scattering of the beam photons inside the cloud. 

A. 2. 3 Approximations 

The OWL detector is located at an altitude of 800 km, this along with the fact 
that the scattering length inside the cloud is on the order of a few km's, implies 
that: Rim w i?om and A± « A±. These variables can then be taken out of the 
integrals and replaced by the approximate values. 

The optical path length due to non-cloud scattering between two points, 
T (nc) ^ re q U j res an integration over the path joining the two points. A significant 
reduction in computation time can be achieved if appropriate approximations 
are used to replace these expressions which appear in the integrals by average 
values which can be taken out of the integrals. In the case of first order scatter- 
ing we have: To™ and . Given the strong forward peak of the scattering 
function we can see that the largest contribution to the integral comes from 
points close to the line joining the source and detector. This allows an approx- 
imation: &i ~ f dss exp(—f3 c s)/ J ds exp(— /? c s) with the integration along the 
line segment from the source to the detector which is contained within the cloud. 
A vector position < ri > can be defined using s%. Now will be replaced 
by an average value and taken out of the integral. 

From the discussion at the beginning of this section we see that aerosols 
and ozone may be ignored in the volume of the cloud. Hence, Tq" w Tq* , 
the optical depth due to Rayleigh scattering. The latter is given by J ds(3a{h) 
where ds is along the line joining rg and t\ and h is the altitude along this line. 
In most cases of interest, the integral can be approximated by |fi — fo| X /3r(/i), 
where h is the height of the midpoint between the two positions. This is due 
to the fact that the cloud thickness is on the order of one to a few km, less 
than the atmosphere scale height of ~ 7 km so (3r does not change much along 
the integration path. This approximation was verified to be accurate to within 
1-5% for a large number of test cases. 

The strong wavelength dependence of Rayleigh scattering implies that the 
calculation should be repeated for each wavelength of interest. However, after 
considering a number of cloud configuration we saw that the result changes by 
less than 20% for wavelengths in the range of 337 nm, and 391 nm. We concluded 
that it would be a reasonable approximation to perform the calculation at a 
wavelength of A = 357 nm, and use the result as an average to be taken out of 
the sum over wavelengths. 

Finally, after the calculation outlines in section A. 2. 2 was implemented for 
the shower's Cerenkov beam, it became apparent that a simple alternative calcu- 
lation which accounts for most of the additional signal received by the detector 
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can be used instead. The strong forward peak of the cloud phase function re- 
sults in that more than 52% of the photons scattered out of a beam are scattered 
forward in a cone of half-angle of 2°. By not subtracting these photons from 
the beam, we in effect calculate the second order scattering of these photons 
at a later stage along the beam propagation when we evaluate the first order 
scattering from the beam at the later stage. The other 48% photons neglected 
in this approximation will have a lesser effect on the detector signal once one 
considers finite time window for the detector pixels. 

A. 3 Simple Cloud Monte Carlo 

A simple method which works well and serves our needs is the Monte Carlo (MC) 
method. MC calculations are valid for clouds of all optical depths, however we 
only employ them for optically thick clouds. Currently our implementation only 
allows for clouds scattering but it can be easily extended to include Rayleigh 
and aerosols scattering. The calculation involves the following steps; 

1. Select photon initial position inside the cloud or at a cloud boundary. 
Also, select the photon direction and a time offset relative to some to. 
The photon direction can be random, for isotropic distribution, or fixed 
in case of a beam. 

2. Propagate the photon by a random step (distance) chosen from an expo- 
nential distribution: exp(— f3 c l), with (3 C the extinction coefficient of the 
cloud 

3. Check if new photon position is inside the cloud. If not then done, if it is 
then continue. 

4. Select a random scattering angle drawn from a distribution which follows 
the clouds phase function. Select a uniform azimuthal angle. Set new 
photon direction. 

5. Goto step 2 

A simple "cloud MC" was developed around this algorithm to calculate the beam 
reflection from the top of a cloud. For a cloud with given cloud parameters a 
large number of photons impinging on the cloud top at a fixed angle (repre- 
senting the beam's zenith angle) is followed through the cloud. As the photons 
emerge from the cloud, either the cloud top (reflected) or cloud bottom (trans- 
mitted) they are added to a set of histograms which record the distributions 
of the locations and time delays of the photons. At the end of the run the 
histograms are saved to file and can be later used by the detector MC. 
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